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Spin systems with frustration and disorder are notoriously difficult to study both analytically and numerically. 
While the simulation of ferromagnetic statistical mechanical models benefits greatly from cluster algorithms, 
these accelerated dynamics methods remain elusive for generic spin-glass-like systems. Here we present a 
cluster algorithm for Ising spin glasses that works in any space dimension and speeds up thermalization by at 
least one order of magnitude at temperatures where thermalization is typically difficult. Our isoenergetic cluster 
moves are based on the Houdayer cluster algorithm for two-dimensional spin glasses and lead to a speedup 
over conventional state-of-the-art methods that increases with the system size. We illustrate the benefits of the 
isoenergetic cluster moves in two and three space dimensions, as well as the nonplanar chimera topology found 
in the D-Wave Inc. quantum annealing machine. 
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A plethora of problems across disciplines map onto spin- 
glass-like Hamiltonians HI]. Despite decades of intense ana¬ 
lytical and numerical scrutiny, a deep understanding of these 
paradigmatic models of disordered systems remains elusive. 
Given the inherent difficulties of studying these Hamiltonians 
analytically beyond mean-field theory as well as the continu¬ 
ous increase of computer power, progress in this field has ben¬ 
efited noticeably from numerical studies. The development 
of efficient Monte Carlo methods such as parallel tempering 
||3] and population annealing 13] has helped in understanding 
these systems at a much deeper level; however, most numeri¬ 
cal studies are still plagued by corrections to finite-size scaling 
due to the small system sizes currently available |0] ■ 

In contrast, simulations of spin Hamiltonians without dis¬ 
order and frustration are comparably simple: Ferromagnetic 
systems have grea^ benefited from the development of clus¬ 
ter algorithms Jj, 13 that help in overcoming critical slowing 
down close to phase transitions. Therefore, the holy grail 
of spin-glass simulations is to introduce accelerated cluster 
dynamics that improve upon the benefits of efficient simula¬ 
tion methods such as population annealing or parallel tem¬ 
pering Monte Carlo. In 2001 Houdayer introduced a semi¬ 
nal rejection-free cluster algorithm tailored to work for two- 
dimensional Ising spin glasses 13] ■ The method updates large 
patches of spins at once, therefore effectively randomizing the 
configurations and efficiently overcoming large barriers in the 
free-energy landscape. Furthermore, the energy of the system 
remains unchanged when performing a cluster move. This 
means that the numerical overhead is very small because the 
rejection rate is zero and there is no need to, for example, 
compute any random numbers for a cluster update. The use 
of these cluster moves made it possible to obtain a speedup 
of several orders of magnitude in two-dimensional systems, 
therefore allowing us to simulate considerably larger system 
sizes. 

While cluster algorithms such as the Swendsen-Wang and 
Wolff ones 13 0 work well for ferromagnetic systems in any 


space dimension because the clusters reflect the spin correla¬ 
tions in the system, this is not the case for algorithms that build 
clusters like the Houdayer cluster algorithm. In this case, the 
clusters do not reflect overlap correlations (313 and cluster 
updates only have an accelerating effect on the dynamics if 
the clusters do not span the entire system or if the comprise 
single spins. This is the case either when temperatures are 
close to zero (small clusters) or when the underlying geome¬ 
try of the problem has a percolation threshold below 50%— 
as is the case in three space dimensions. Updating such a 
system-spanning cluster amounts to swapping out both repli¬ 
cas, thereby not randomizing the configurations. This means 
that while the method works in principle, it does not really 
provide any simulational benefit. As such, Houdayer cluster 
moves work, in principle, only for models where the percola¬ 
tion threshold is above 50%, as is the case in two-dimensional 
Ising spin-glass Hamiltonians. One way to remedy this situa¬ 
tion is to increase the percolation threshold artificially, e.g., by 
diluting the lattice llOtl . However, this is often not desirable 
and is highly dependent on the problem to be studied. 

Here, we show that Houdayer-like cluster moves can be 
applied to spin systems on topologies where the percolation 
threshold is below 50%, provided that the interplay of tem¬ 
perature and frustration prevents clusters from spanning the 
whole system. We therefore introduce isoenergetic cluster 
moves for spin-glass-like Hamiltonians in any space dimen¬ 
sion. These rejection-free cluster moves accelerate thermal¬ 
ization by several orders of magnitude even for systems with 
space dimensions larger than 2. We show that the inherent 
frustration present in spin-glass Hamiltonians prevents clus¬ 
ters from spanning the whole system for temperatures below 
the characteristic energy scale of the problem. As such, spin- 
glass simulations can be sped up considerably in the hard-to- 
reach low-temperature regime of interest in many numerical 
studies. 

The fact that the isoenergetic cluster moves are rejection 
free and leave the energy of the system unchanged is also 
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of great importance to any heuristic based on Monte Carlo 
updates to compute ground-state configurations of spin-glass- 
like Hamiltonians. For example, the convergence of simulated 
annealing ifTlI] can be considerably improved by adding isoen- 
ergetic cluster moves at each temperature step. Because the 
moves change the spin configurations but leave the energy of 
the system intact, the approach has the potential to “tunnel” 
through energy barriers, thus improving overall convergence. 

We first introduce the benchmark model, followed by a 
short description of the Houdayer cluster algorithm and an 
outline of our isoenergetic cluster algorithm. Results in two 
and three space dimensions, as well as on the nonplanar 
chimera topology lll2ll are presented. 

Benchmark model and observables. —The Hamiltonian of 
a generic Ising spin glass is defined by "H = JijSiSj, 

where Si € {±1} represent Ising spins and N is the total num¬ 
ber of spins. In this study the interactions Jij are selected from 
a Gaussian distribution with mean zero and variance = 1. 
Because we are only interested in highlighting the improved 
thermalization by adding isoenergetic cluster moves, we mea¬ 
sure the average energy per spin defined via \{'H)\/N, as well 
as the link overlap = (l/.^Vb) . Here, 

(• • •) represents a Monte Carlo average, the superscripts rep¬ 
resent two replicas of the system, [• • • ] indicates an average 
over the disorder, and iVb is the number of bonds in the sys¬ 
tem. Using Gaussian disorder, one can equate the internal 
energy per spin to the internal energy computed from the link 
overlap S, i.e.. 




( 1 ) 


To test that the system is thermalized, we thus study the time- 
dependent behavior of 
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sites are flipped with probability 1, irrespective of their ori¬ 
entation. The method can be implemented in a very efficient 
way because sites are added to the cluster with probability 1 
and the cluster updates are rejection free. To ensure ergod- 
icity, the cluster move is combined with standard single-spin 
Monte Carlo updates. Summarizing, one simulation step us¬ 
ing the HCA consists of the following steps: 

1. Perform one Monte Carlo sweep (TV Metropolis up¬ 
dates) in each replica. 

2. Perform one Houdayer cluster move. 

3. Perform one parallel tempering update for a pair of 
neighboring temperatures. 

Note that the last step is not necessary; however, the combi¬ 
nation of the HCA moves and parallel tempering (PT) updates 
improves thermalization considerably and represents the stan¬ 
dard modus operandi. 

In theory, the efficiency of the HCA depends strongly on the 
percolation threshold of the desired topology to be simulated. 
Because spins are added to the cluster with probability 1, if 
the percolation threshold of the studied lattice is below 50%, 
then the cluster might span the entire system and an update 
will not yield a new configuration. This is the reason wlw 
the HCA is claimed to only work in two space dimensions ITl] 
where the percolation threshold is above 50% (see also Fig.[T] 
top panel). 

Isoenergetic cluster algorithm. —Our proposed isoenergetic 
cluster moves are closely related to the HCA. We begin by 
simulating two replicas with the same disorder at multiple 
temperatures. The cluster moves alone are not ergodic, so, 
again, these must be combined with simple Monte Carlo up¬ 
dates. One simulation step using isoenergetic cluster moves 
consists of the following steps: 


When A —0, the bulk of the disorder instances are thermal¬ 
ized il4ll . Simulation parameters are listed in Table U 

Reminder: Houdayer cluster algorithm. —The Houdayer 
cluster algorithm (HCA) ||3l is an efficient algorithm to study 
two-dimensional Ising spin glasses at low temperatures where 
thermalization is slow. It is similar to replica Monte Carlo 
1 151], but with the difference that both replicas are at the same 
temperature. By allowing large cluster rearrangements of con¬ 
figurations, the HCA improves thermalization by efficiently 
tunneling through configuration space. 

The algorithm works as follows: In the HCA, two indepen¬ 
dent spin configurations (replicas) are simulated at the same 
temperature. The site overlap between replicas (1) and (2), 


Qi = is calculated. This creates two domains in q 

space: sites with qi = 1 and qt = —1. Clusters are defined 
as the connected parts of these domains in q space. One then 
randomly chooses one site with qi = —1 and builds the clus¬ 
ter by adding all of the connected spins in the domain with 
probability 1. When no more spins can be added to the cluster 
in q space, the spins in both replicas that correspond to cluster 


1. Perform one Monte Carlo sweep (TV Metropolis up¬ 
dates) in each replica. 

2a. If the number of cluster sites with qi = —1 is greater 
than TV/2, then all the spins in one of the configura¬ 
tions can be flipped (because of spin-reversal symme¬ 
try), thus reducing the cluster size while leaving the en¬ 
ergy unchanged. 

2b. Perform one Houdayer cluster move for all tempera¬ 
tures T < J. 

3. Perform one parallel tempering update for a pair of 
neighboring temperatures. 

The main difference thus lies in applying cluster moves to 
a carefully selected set of temperatures where the isoenergetic 
cluster moves (ICMs) are efficient (steps 2a and 2b) because 
clusters do not percolate, as well as reducing cluster sizes 
and thus the numerical overhead by exploiting spin-reversal 
symmetry (step 2a) SB. For example, in the case of the 
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chimera lattice the overhead of the ICM over PT is approxi¬ 
mately 25% and is roughly independent of the system size for 
the studied N. However, the overhead for the HCA over PT is 
at least 50% and grows with increasing system size. 

Figure [T] shows the fraction of spins with negative overlap 
(i.e., the fraction of potential cluster sites) as a function of 
temperature T for different system sizes N and on three differ¬ 
ent topologies. The top panel of Fig.[T]shows data in two space 
dimensions where the percolation threshold is pc ~ 0.592 
1 181] (the solid horizontal line). As such, for all temperatures 
simulated, the fraction of cluster sites is below the percola¬ 
tion threshold and saturates at 50% for T —?> oo. This means 
that isoenergetic cluster updates are efficient for all tempera¬ 
tures studied because the clusters never percolate. Naively, 
one would expect that in higher space dimensions clusters 
percolate for all T’s. This is, however, not the case due to 
the frustration present in spin glasses, as can be seen for the 
chimera topology (the center panel of Fig. [U or in three space 
dimensions (the bottom panel of Fig. [T]). For increasing sys¬ 
tem size the fraction of cluster sites converges to a limiting 
curve that crosses the percolation threshold (the horizontal 
solid lines) at approximately T k, J = \. This means that, 
for all T > J, clusters percolate and the cluster updates are 
just numerical overhead without any advantage to the simu¬ 
lation. However, for T < J the fraction of cluster sites lies 
below the percolation threshold. This means that perform¬ 
ing cluster moves in this temperature regime should improve 
thermalization. Note that it is a coincidental property that for 
three-dimensional Ising spin glasses Tc 1 |l9ll. i.e., that 
cluster moves can be applied to any T <Tc 1201]. 

When the interactions Ja are drawn from a Gaussian dis¬ 


tribution, the ground state is unique. As can be seen in Fig.[Tl 
the fraction p of spins potentially in a cluster also approaches 
zero for T —?> 0; i.e., both replicas are in the ground state for 
low enough T. Therefore, the cluster is composed of no sites 
or the entire lattice. In the case of disorder distributions that 
yield a highly degenerate ground state, such as is the case for 
bimodal disorder, it is possible to continue to have clusters at 
zero temperature. It is thus possible to efficiently hop around 
the ground-state manifold by applying cluster moves to low- 
lying or even zero-temperature states, although this might not 
be ergodic. We do emphasize, however, that if clusters are too 
small, then the isoenergetic cluster moves also become inef¬ 
fective. Therefore, plotting p as was done in Fig.[T]is essential 
in determining the efficiency and applicability of the method. 

Benchmarking results .—Figure |2] shows A [Eq. ©j as a 
function of Monte Carlo time (measured in lattice sweeps) 
t = 2^. The top panel of Fig. |2] shows data in two space 
dimensions for simulations using isoenergetic cluster moves 
(PTh-ICM) and vanilla PT Monte Carlo for N = 1024 spins 
at T = 0.212. Once A ~ 0, we deem the system thermal- 
ized. Clearly, the inclusion of cluster moves—as can also be 
expected from the results of Houdayer—show an improved 
thermalization. The center panel of Fig. |2] shows data on the 
chimera topology with N = 1152 spins and T = 0.212, 
where the HCA is not expected to show any improvement 
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FIG. 1: (color online). (Top panel) Fraction of spins p of potential 
cluster sites as a function of temperature T for different system sizes 
N in two space dimensions (2D). The horizontal line represents the 
percolation threshold of a two-dimensional square lattice, i.e., pc ~ 
0.592 fl^ . Because p —>■ 0.5 for T — l oo, for all T clusters do 
not percolate, which is why the HCA is efficient in two-dimensional 
planar geometries. (Center panel) p as a function of temperature T 
for different system sizes N on the chimera topology. The horizontal 
line represents the percolation threshold of the nonplanar chimera 
topology, namely Pe « 0.387 computed here using the approach 
developed in Ref. (see the Supplemental Material). For T > J = 
1 clusters percolate and cluster updates provide no gain. (Bottom 
panel) p?&di function of temperature T for different system sizes N 
in three space dimensions (3D). The horizontal line represents the 
percolation threshold of the three-dimensional cubic lattice (pc « 
0.311 1^ ). For T > J = 1 clusters percolate. In all panels, error 
bars are computed via a jackknife analysis over configurations and 
are smaller than the symbols. 
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FIG. 2: (color online). (Top panel) A [Eq. m as a function of 
simulation time t — 2^ measured in Monte Carlo sweeps in two 
space dimensions (2D) for N = 1024 and T — 0.212. Simula¬ 
tions using vanilla PT thermalize at at least 2^® Monte Carlo sweeps, 
whereas with the addition of ICMs thermalization is reduced to ap¬ 
proximately 2^® Monte Carlo sweeps. This means approximately 2 
orders of magnitude improvement. (Center panel) A as a function 
of simulation time f = 2® measured in Monte Carlo sweeps for an 
Ising spin glass on chimera with N = 1152 spins at T = 0.212. 
Simulations using PT thermalize at approximately 2^® Monte Carlo 
sweeps, whereas the addition of ICMs reduces thermalization to 2^® 
Monte Carlo sweeps. Again, approximately 2 orders of magnitude 
speedup. (Bottom panel) A as a function of simulation time t = 2® 
measured in Monte Carlo sweeps in three space dimensions (3D) for 
N = 1728 and T = 0.42 ~ 0.43Tc. Using standard PT, the sys¬ 
tem thermalizes approximately after 2^® Monte Carlo sweeps. This 
time is reduced to ~ 2®° Monte Carlo sweeps when ICMs are added. 
In all panels, error bars are computed via a jackknife analysis over 
configurations. 


TABLE I: Parameters of the simulation in two space dimensions 
(2D), three space dimensions (3D), and on the chimera (Ch) topol¬ 
ogy. Eor each topology simulated and system sizes N, we compute 
Nss. disorder instances and measure over 2® Monte Carlo sweeps 
(and isoenergetic cluster moves) for each of the 2Nt replicas. Tmin 
[Tmax] is the lowest [highest] temperature simulated, and Nt is the 
total number of temperatures used in the parallel tempering Monte 
Carlo method. Isoenergetic cluster moves only occur for the lowest 
Ac temperatures simulated (determined from Pig.[T}- 

N Nsa b Tmin Tmax Nt Nc 

2D 256,576, 1024 W 22 0.2120 1.6325 30 ^ 

Ch 128, 288, 512, 800, 1152 10'‘ 22 0.2120 1.6325 30 19 

3D 64, 216, 512, 1000, 1728 1.5 10'‘ 23 0.4200 1.8000 26 13 


over PT due to pc < 0.5. As can be seen, our ICM clearly 
improve thermalization in comparison to PT by at least 2 or¬ 
ders of magnitude, an amount that grows with increasing sys¬ 
tem size. Finally, the bottom panel of Fig. |2] shows A as a 
function of simulation time in three space dimensions with 
N = 1728 spins and T = 0.42 <C T].. Although not as im¬ 
pressive as with the chimera topology, we see a speedup of 
approximately one order of magnitude—an amount that again 
grows with increasing system size. 

Finally, Fig. [3 shows the ratio of the thermalization time 
using PT and using PTi-ICM for different topologies at the 
lowest simulation temperature (see Table HI as a function of 
the system size N. In all cases, the speedup increases with 
increasing system size, therefore illustrating that the addition 
of isoenergetic cluster moves greatly improves thermalization. 
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FIG. 3: (color online). Ratio between the approximate average ther¬ 
malization time of PT and PT-l-ICM for different topologies at the 
lowest simulation temperature (see Table |I| as a function of system 
size N. In all cases the speedup increases with increasing system 
size. Note that thermalization times have been determined by eye. 

Summary .—We have presented a rejection-free cluster algo¬ 
rithm for spin glasses in any space dimension that greatly im¬ 
proves thermalization. By restricting Houdayer cluster moves 
to temperatures where cluster percolation is hampered by the 
interplay of frustration and temperature, we are able to ex¬ 
tend the Houdayer cluster algorithm for two-dimensional spin 
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glasses to any topology or space dimension. Our standard im¬ 
plementation of the cluster updates represents only a minor 
overhead compared to the thermalization time speedup 
obtained from the isoenergetic cluster moves—a speedup that 
increases with the system size UM . 
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Supplemental Material: Percolation threshold of the 
chimera lattice. —In percolation theory, the relative size of 
the largest cluster of a lattice (smax) plays the role of 
an order parameter. Above the percolation threshold pc, 
(sinax) —1, whereas below pc it approaches 0. One can 
define a dimensionless Binder ratio 12ill via g = (1/2) [3 — 
(('Smax)/('Smax)^)] determine the percolation threshold p^ 
to high precision. Here, (• • •) represents an average over 
randomly-generated configurations with a site probability p 
on a chimera lattice with N sites. Close to criticality, g ~ 
l/pc)], i.e., whenp = Pc data for different 
system sizes cross. Figure |4] shows a finite-size scaling plot 
of the data for g. We estimate for the percolation threshold of 
the chimera lattice pc = 0.3866(3). 
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